How urban form impacts flooding

Urbanization and climate change are contributing to severe flooding globally, damaging infrastructure, disrupting economies, and undermining human well-being. Approaches to make cities more resilient to floods are emerging, notably with the design of flood-resilient structures, but relatively little is known about the role of urban form and its complexity in the concentration of flooding. We leverage statistical mechanics to reduce the complexity of urban flooding and develop a mean-flow theory that relates flood hazards to urban form characterized by the ground slope, urban porosity, and the Mermin order parameter which measures symmetry in building arrangements. The mean-flow theory presents a dimensionless flood depth that scales linearly with the urban porosity and the order parameter, with different scaling for disordered square- and hexagon-like forms. A universal scaling is obtained by introducing an effective mean chord length representative of the unobstructed downslope travel distance for flood water, yielding an analytical model for neighborhood-scale flood hazards globally. The proposed mean-flow theory is applied to probe city-to-city variations in flood hazards, and shows promising results linking recorded flood losses to urban form and observed rainfall extremes.

properties of streams 14 , Q = AR 2=3 ffiffiffi α p =n M , where flow rate Q is related to Manning's roughness n M , flow cross-section and hydraulic radius A and R, and the bottom slope α, respectively.The Manning equation and other empirical models (e.g., Chezy) have proven instrumental and robust for designing flood channels, and here we seek an equivalent for urban flood plains.
Our inspiration for an urban flooding mean-flow theory is taken from how complex systems are approached in statistical mechanics.Urban forms epitomize a complex system much akin to granular media 15 , disordered porous solids 16 , glassy systems 17 , and complex fluids 18 , to name only a few.When viewed deterministically, flooding through a complex urban form depends on the location of each obstruction.However, in an average sense, the flow behavior should be governed by a few characteristic state variables reflective of the urban form, e.g., building density, and indicators of how clusters of buildings are configured.This averaging is motivated by Boltzmann's treatise that proves the equilibrium properties of non-interacting classical particles remain independent of their location and instead rely on the system's volume-a fundamental thermodynamic state variable 19 .To reduce the complexity of urban flooding, we need a large ensemble of refined flow observations to build our mean-flow theory.Experimental studies 20 and field measurements 21 have proven insightful when studying the role of obstructions on urban flooding, yet to date, they remain strictly limited to extremes of simplified obstruction arrangements or case-specific urban forms.Hence, computational modeling is the only viable approach to collect statistically significant data on urban flooding dynamics over a representative spectrum of urban forms.
The flow structure and boundary layers of urban flooding can be simulated with Reynolds-averaged Navier-Stokes, large-eddy simulation, and direct numerical simulation 22 .However, in pursuit of computational efficacy, we rely on a shallow-water hydrodynamic model 23 and produce thousands of fine-resolution, building-resolved simulations for different realizations of urban form.Shallow-water theory builds on the assumption of hydrostatic pressure and nearly horizontal flow, and shallow-water models have been validated at neighborhood and larger scales by extensive laboratory and field studies and urban flood modeling 20,[24][25][26][27][28][29] .While urban flooding can involve high levels of complexity, such as flows into basements, subway systems, and usually involves a mix of underground flows through storm pipes and overland flow along streets 30 , we are specifically interested here in the depth and velocity of overland flow through cities caused by the amount of overland flow in response to hazard drivers (e.g., rainfall, overbank flows, stormwater surcharge).Hence, our simulation approach involves measuring the steady-state depth and velocity of flooding for a given flow rate under a wide range of urban forms.
In this work, we assemble an ensemble of urban flood simulations and reduce the complexity of urban flooding through dimensional analysis of simulation results and consideration of conservation laws into a master curve that constitutes a mean-flow theory.This theory establishes a mathematical equation that relates flooding attributes-notably depth and intensity defined as the product of flood depth and velocity-to the ground slope, urban form, and flow rate encountered across cities.Moreover, we seek a model that explains the global variability of urban flooding hazards at neighborhood and larger scales and provides planning-level support for more flood-resilient cities.

Results and discussion
Mean-flow theory for urban flooding City building arrangements are highly varied globally, as indicated in Fig. 1 (Supplementary Note I and Materials and "Methods" Section A).For example, the alignment of buildings in Chicago, USA (Fig. 1c) is structured, while Lagos, Nigeria is unstructured (Fig. 1h).A meanflow theory of urban flooding applicable to global variations in building layouts (e.g., Fig. 1(a-h)) calls for a systematic parameterization of urban form, and here we draw from studies of the physics of phase transition and disordered materials 31,32 , where the Mermin order parameter 33 (χ) has been used to quantify the spatial degree of order/disorder of a system of particles.At the scale of individual particles (or buildings), the Mermin order parameter determines the symmetry of c n nearest neighbors to the reference particle 34 .For instance, for building arrangements with square symmetry, c n = 4, and those with hexagonal symmetry feature c n = 6.Furthermore, the local order parameter can be averaged to extract a collective disorder attribute representative of an urban region, as follows: where θ kjk 0 is the angle between the central building j and its neighbors k and k 0 , where k 0 is a fixed neighbor and k permutes over all possible neighbors of the jth building.When χ c n is greater than 0.9, the synthetic urban form exhibits a pronounced symmetry and resembles a planned urban layout, e.g., Chicago, Fig. 1c.As disorder increases, the symmetry diminishes, and the synthetic urban form features characteristics of an organically nucleated and grown urban environment, e.g., Jakarta, Fig. 1e.Here, we adopt the order parameter concept to quantitatively characterize the disorder of buildings within a city alongside the urban porosity ϕ, which characterizes the amount of pore space and is inversely related to the building density.Ensembles of synthetic city urban forms at a given porosity and disorder are developed using hybrid reverse Monte Carlo 35,36 (HRMC); see Materials and "Methods" Section B and Supplementary Notes II-III for detail.Furthermore, to support shallow-water modeling and specifically the characterization of steady-state, uniform flow conditions, a rectangular urban domain is created with a length (L), width (W), longitudinal slope (α), and a Manning coefficient for bottom resistance (n M ), and the domain is populated with N square buildings of side length (D).Due to HRMC's computational cost, we only generate building patterns in a unit cell of length l u and repeat the unit cell to cover the entire channel length.The chosen level of disorder controls the irregularity of building locations (see Supplementary Fig. S1), and the urban porosity follows as ϕ = 1 − ND 2 /(WL).To characterize flooding, the shallow-water equations are solved on a fine-resolution, unstructured computational mesh (ParBreZo), where mesh edges are aligned with building walls to enforce a no-flow boundary condition 23,28 , see Materials and "Methods" Section C, Supplementary Figs.S2-S3.
Our simulations resolve flooding within the pore space between buildings, and solutions reflect the bulk effect of building form drag.Moreover, using a Godunov-type finite volume scheme for solving the shallow-water equations allows us to resolve supercritical flows, subcritical flows, and flows with hydraulic jumps 23 .Solutions include localized super-critical flows (Fr > 1) along steep slopes, but average flows across neighborhoods are found to be largely sub-critical (Fr < 1).
As has been observed in experiments 37 , simulations show that a hydraulic jump occurs when localized supercritical flows encounter obstructions 37 .
Figure 2a-f show the color map of non-dimensionalized flood velocity in the downslope direction (uD 2 /Q) as a function of flow velocity (u), flow rate (Q), and building size (D), for neighborhoods of varying porosity and disorder.Relative to square symmetry at ϕ = 0.6, we observe the choking of flow with decreasing χ 4 .This indicates that increasing disorder restricts the flood flow in urban forms that possess a near square symmetry (regular layout), but has an opposite effect in synthetic urban models with nearly hexagonal symmetry (or a staggered layout) representative of flow paths along diagonal directions.In particular, at ϕ = 0.6, the disorder opens new flood flow pathways through clusters of buildings much akin to flow pathways through disordered granular media 38 .
These qualitative observations suggest that both the symmetry and irregularity of street grids in cities contribute to the tortuosity and constriction of flow paths at neighborhood and larger scales, which in turn shape the distribution of flood depths and velocities responsible for safety risks and losses 39 .To characterize how the symmetry and irregularity interact with ground slope to control the neighborhoodscale flood levels, we leverage the non-dimensionalized chord length distribution Pðl 0 c ðθÞÞ, which represents the distribution of randomly placed line segments fitting within pore space between two building walls 40 in all directions (θ), see Supplementary Note IV and Supplementary Fig. S4 for detail.For simplicity, we here focus on chord length in the direction of the steepest descent.The chord length represents the distance over which a fluid particle can accelerate under the influence of gravity and bottom shear before being blocked by an obstruction, and Fig. 2(g-l) show the distribution of Pðl 0 c Þ for the building patterns shown in Fig. 2(a-f), respectively.The floodplain with perfect square symmetry exhibits a bimodal distribution of long and short chord lengths corresponding to the unit cell's length and nearest neighbor distance.By introducing disorder, additional peaks appear at intermediate values of chord length.The perfectly hexagonal floodplain also displays a bimodal distribution and behaves similarly; however, the large chord length corresponds to the second nearest neighbor.A representative chord length for any urban form follows with the expectation of l 0 c excluding the first and second nearest neighbor chords that do not effectively contribute to the flood drainage due to stagnation effects.Focusing on chord lengths beyond three building sizes, we define the normalized average effective chord length as: where l 0 c = l c =D and l u is the length of the unit cell.The corresponding l c and average non-dimensionalized velocity for each distribution also appear in Fig. 2g-l.We observe qualitatively that the average flood velocity is directly correlated with l c , meaning that long l c results in higher average velocities regardless of the underlying urban symmetry.
To reduce the observed complexity in urban form-induced flood behavior, we consider the mass and momentum conservation laws for dense urban forms, where the bottom shear is negligible compared to the distributed form drag.In a spatially large enough control volume that sufficiently represents the flow heterogeneity, the average gravitational force ( F g = g hAϕ cosðαÞsinðαÞ) balances the resultant form drag ( F D = C D A p ujV j) in the steady-state limit, Fig. 3(a).In this case, while gravity (g) drives the flow with average vertical height ( h) through available urban pore space (Aϕ), the drag resists the flow.This resistive force is proportional to the drag coefficient (C D ), the so-called projected area (A p = hW p cosðαÞ) 41 , and spatially averaged velocity The average perpendicular component of flow velocity (v) is zero, and its magnitude should be smaller than u, meaning ujV j≈u 2 .Invoking mass conservation (q = uhϕ cosðαÞ), we can eliminate the drag dependence on velocity and estimate this force using average height, hϕ 2 cosðαÞÞ.Therefore, through force balance ( we can estimate the average flood height as follows, The main challenge in determining the average flood height pertains to the exact measurement of projected width (W p ), an extensive quantity that depends on urban form.Therefore, Aϕ 3 /W p can be viewed as a characteristic length scale whose size is proportional to the building size, Aϕ 3 /W p = D/Π 2 , where similar to C D is a dimensionless function of the building arrangements (r i ) in the urban layout.For the maximum 10% slope considered in this study, the approximation of cosðαÞ ffiffiffiffiffiffiffiffiffiffiffiffiffi sinðαÞ p / ffiffiffi α p introduces less than 1% error.
The average flood height in Eq. ( 3) is deterministic if the locations of all buildings are known.However, when viewed through statistical mechanics, the flood height should, in principle, be only dependent on a small number of statistical attributes of the floodplain rather than the unique locations of all buildings.While this is an unconventional concept in the context of urban flooding, it is a well-established concept for studying the thermodynamics of granular media and glassy systems 17 .In this sense, the ensemble-averaged dimensionless flood height can be thought of as a hydrodynamic state function of the city that depends on hydrodynamic state variables ϕ, χ, and l c , which yields the following dimensionless relation: where 〈.〉 denotes the ensemble average and modified drag coefficient Interestingly, we derive the same relation through considerations of dimensional homogeneity, see Supplementary Note V. Regardless of the derivation method, the theory suggests the ensemble-averaged dimensionless flood height (h hi ffiffiffiffiffiffi ffi gD p =q) should depend on dimensionless geometric factors, the bottom slope, and the drag coefficient that depends on the porosity, order parameter, and average chord length in the flow direction.To test this theory, we perform one-factor-at-a-time flow simulations varying slope and flow rate per width for different urban forms with various porosity and order parameters, Fig. 3(b-c).Corroborating with the proposed theory, the numerical simulations show h hi scales linearly with q and is inversely proportional to the square root of the bottom slope, 1= ffiffiffi α p .We note that while our theory and Manning equation show h hi / 1= ffiffiffi α p , they differ in predicting the exponent in h hi / q m scaling relation.In an empty wide rectangular channel, where flow is driven by gravity and bottom shear, the Manning equation predicts h hi / q 3=5 .However, as discussed before, in a densely-packed urban flood channel, the form drag balances the gravity leading to a new scaling relation, i.e., h hi / q 1 .Therefore, the proposed change in exponent m is intimately related to the change in underlying physical processes governing momentum transfer.
To further verify the proposed theoretical framework, we need to demonstrate that C 0 D ðϕ,χ, l c Þ can solely explain the distributions of flood depth over a range of flow simulations with varying 0.4≤ϕ≤1, 0.3≤χ 4 , χ 6 ≤1, α, and q.In an attempt to collapse the results, we first neglect the impact of urban order, i.e.,C 0 D = f 1 ðϕÞ / ð1 + a c n ϕÞ s c n =ϕ p c n , where a c n , s c n and q c n are specific to square and hexagonal symmetries as denoted by the coordination number c n .This approach reduces the search for the state function to an optimization problem.Regardless of the underlying symmetry, s c n ≈1, p c n ≈0 and a c n converges to -1, meaning f 1 ≈ (1 − ϕ), better known as the packing fraction in the physics of granular media.This parameter mirrors the urban building coverage, which is suggested to be the most influential factor affecting flood depth via a sensitivity analysis 12 .The coverage parameter also appears in flow through channels covered by vegetation and is shown to scale flow resistance 42 .However, the packing fraction collapses data only at high porosity and becomes increasingly insufficient in dense textures, where the disorder becomes a controlling factor, see Supplementary Fig. S8a.
To account for the disorder effects at low porosity, we assume C 0 D = f 1 ðϕÞ × f 2 ðχÞ.For simplicity, we focus on linear functions through the Taylor expansion, C 0 D ≈c c n ð1 À ϕÞð1 + b c n χ c n Þ, where b c n and c c n are fitting parameters.In particular, we arrive at the following relation: where slopes c 4 and c 6 are 5.31 and 2.54, and b 4 and b 6 are -0.5 and 1.1, respectively, as noted in Fig. 4a.Interestingly, we note c 4 /c 6 ≈ 2 and b 4 / b 6 ≈ − 1/2.As qualitatively observed in Fig. 2, the introduction of disorder into square and hexagonal urban forms has an opposite effect on flood behavior.This is now quantitatively confirmed by the change in the b c n sign.We also note the proposed linear relationships do not capture the flood depth for floodplains without obstructions.The relevant dimensionless quantity in the limit of high ϕ is h hi 5=3 =ðqn M Þ, which is confirmed by the convergence of simulation data when ð1 À ϕÞð1 + b c n χ c n Þ= ffiffiffi α p tends to zero, see Supplementary Fig. S8b.
We can now examine the impact of urban form on flood depth in synthetic city neighborhoods.Leveraging Eq. ( 5) at constant porosity, the flood depth ratio between square-and hexagon-like patterns becomes c 4 (1 + b 4 χ 4 )/(c 6 (1 + b 6 χ 6 )).This ratio is ≈ 1/2 for perfectly ordered cities, meaning flood depth is twice as severe in hexagonal cities.As the disorder increases to the point that the underlying patterns become indistinguishable, i.e.,χ 4 = χ 6 = 0.5, the ratio becomes ≈ 1, and as expected, there is no difference in flood depth.In this scenario, the disorder increases flood depth by ≈ 1/2 in a square urban pattern, while it decreases the flood depth by ≈ 1/4 in hexagonal form.
The inherent drawback of Eq. ( 5) as a mean-flow theory is that it does not provide a universal relation between dimensionless quantities due to dependence on the underlying symmetry.As noted in Fig. 2, the flood depth is associated with open flow pathways characterized by l c .This prompts proposing C 0 D ðϕ, Focusing on an f 3 function of the form l c β , we find through optimization that β ≈ − 1/2.As shown in Fig. 4b, this heuristic approach collapses the bifurcated results into a linear relationship as follows: where d 0 = 0.07 and d 1 = 5.62.Eqs. ( 5)-( 6) are valid in the range of 0.4≤ϕ considered here.Eq. ( 5) becomes increasingly inaccurate at lower porosity as flood depth should increase non-asymptotically at ϕ = 0. Eq. ( 6), however, becomes singular at ϕ = 0, since l c tends to zero in this regime.This equation constitutes a mean-flow theory of urban flooding that relates the average flood depth to urban form attributes.Note that increasing slopes and increasing chord length act to decrease flood depth due to increases in flood velocity.Hence, the mean-flow theory predicts shallow, high-velocity flows along streets oriented in the downslope direction, an especially dangerous flow regime that has been described as ultrahazardous flooding 43,44 .

Global Flood Hazard Distribution based on the Mean-Flow Theory
The proposed mean-flow theory is now applied to examine the city-tocity variability of urban flooding hazards; see Materials and "Methods" Section D. Selecting twenty cities across five continents, we analyze urban porosity and chord length in the direction of the steepest slope at a km 2 grid resolution as shown in Fig. 5(a-b).Similar to the synthetic floodplains, we observe a correlation between urban porosity and normalized average chord length, Fig. 5c.The denser cities feature smaller chord lengths, as evident in San Francisco and Lagos.However, a notable outlier to this trend is São Paulo, where long streets in densely developed neighborhoods create long corridors for flood passage.Loosely-packed cities such as Virginia Beach also exhibit longer chord lengths due to the sparsity of flow obstructions.
To characterize cities in terms of their flood hazards, we define two dimensionless flood hazard indicators: flood depth and intensity; see Supplementary Note VI.For the dimensionless flood depth hazard H * , we use our mean-flow theory in Eq. ( 6), which is independent of the underlying symmetry and more applicable to real-world complex urban forms.Flow rates required by the mean-flow theory, q, are estimated as the product of an extreme event precipitation intensity and a length scale for runoff accumulation, taken here to be 1 km.We then divide h hi by h ref , the reference flood depth in an obstruction-free channel with α ref = 0.1 and q ref = 0.001 m 2 /s, which are respectively representative of average slopes and precipitation rates observed in our twenty cities.We also define the dimensionless flood intensity as the average product of flood depth and velocity normalized by the flood intensity in an obstruction-free channel inundated by q ref , P * = hhui= ðhuÞ ref .This dimensionless quantity is controlled by the inflow rate at the steady state, meaning P * = q/(q ref ϕ).This indicates the normalized flood intensity at the urban scale is independent of its form and bottom slope, and controlled by porosity and flow rate only; see Supplementary Fig. S9.
How susceptible are cities to urban flood hazards based on urban form alone?We can now address this question by computing the dimensionless flood depth and intensity given the urban porosity, chord length, and slope for each city, as shown in Fig. 5d.Note that with q computed with a 1 km runoff accumulation length scale for all cities, differences in flood depth and intensity are completely attributable to differences in urban form and precipitation intensity, and not to differences arising from larger scale drainage patterns or urban drainage infrastructure that contribute to the accumulation of urban flood flows.Most cities are clustered toward the low-hazard region, i.e., low flood depth and intensity.Chicago's urban form is conducive to high flood depth due to its relatively low porosity and flat terrain.San Francisco faces the opposite hazard, where the flood intensity overshadows flood depth due to the steep topography and low urban porosity.High-hazard cities are marked by high levels of flood depth and intensity and include Tokyo, Lagos, Jakarta, Karachi, and São Paulo.Compared to other cities in this study, these cities exhibit a very low porosity and mild slope.In the case of São Paulo, there are also long and steep flow pathways characterized by high chord length and bottom slope.
We reviewed the Emergency Events Database (EM-DAT) 45 for extreme events that occurred within the 20 cities considered in this study and used archives of historical precipitation and a simple rainfallrunoff approximation to estimate the inflow rate, q, see Materials and "Methods" Section D. Figure 5(e) maps hazard severity for sixteen extreme events classified as flash floods using our dimensionless flood depth and intensity.São Paulo, Jakarta, Karachi, Toyko, and Houston emerge as the cities that have experienced the most intense flash flood hazards, based on dimensionless flood depth, intensity, or a combination of the two.Moreover, the occurrence of multiple events in São Paulo, Jakarta, and Karachi suggests that these cities are amongst the most hazardous cities in the world for flash floods based on a combination of extreme event frequency and urban form.
Validating the proposed mean-flow theory for describing city-tocity variations in urban flooding is extremely difficult because extreme event flooding of street grids is not systematically documented with ground-based sensors and satellite-based remote sensing methods for flood inundation do not perform well in urban areas 46 .EM-DAT includes metrics of economic losses, which we expect to be associated with precipitation, flood depths, and flood velocity.However, flood losses depend on many physical and social vulnerabilities in addition to hazard severity.Thus any attempt to link hazard severity to economic impacts is expected to be highly uncertain.We focused on five flash flooding events from EM-DAT with economic loss data within the twenty selected cities of this study.Due to the economic variability of the impacted cities, we normalize the reported damages by the respective country's gross domestic product and re-scale the damage values between zero and unity.Figure 5(f) reveals a correlation between the normalized monetary damages, L, and a linear combination of flood hazard indicators as follows: where This relationship provides a rather intuitive yet insightful finding that flood damages are linearly related to q and an urban form factor, F .The positive correlation between L and F is in general alignment with established theories and models for flood damage, including widely used depth-damage curves [e.g., 47 ].If we set κ = 0 in our damage model, Eq. ( 8), it is simplified to a linear depth-damage curve.However, we find that the reported EM-DAT damages cannot be solely explained by flood height, which necessitates the consideration of flood intensity as a damage-driving factor.The Houston (2016) event appearing in Fig. 5f was the most costly of the flash flood events recorded across the cities examined for this study.However, the most costly extreme events in recorded history were all classified by EM-DAT as tropical cyclones, not flash floods, and the top six events as of 2019 all occurred in the U.S. 2 .When we compared tropical cyclone losses and our estimates of urban flood hazards from rainfall, we did not find a correlation.Indeed, damage from tropical storms results from several hazard drivers, including wind, storm surge, and rainfall 2 .
We acknowledge several limitations of our approach, including the representation of urban forms with arrays of square obstructions as opposed to more organic forms, assumed uniformity in the runoff accumulation length scale, the approximation of urban flooding as a two-dimensional flow without fully resolving the turbulent boundary layers that contribute to form drag, and the limited number of cities sampled for the development of our mean-flow theory.Additional research could help with understanding the importance of each of these methodological factors relative to the development of a meanflow theory for urban flooding.Nevertheless, given the range of urban porosity, order values, and floodplain slopes considered herein, our approach systematically considers an immense spectrum of urban forms.Indeed, we successfully reduce the complexity of urban flood flow data into a master curve that provides an intuitive explanation of the distribution of flood hazards globally across cities, and that aligns with observations of flood losses.
The mean-flow theory presented here thus offers a foundational approach for examining how urban form affects the physical vulnerability of cities to flooding.Application of the mean-flow theory at neighborhood and larger scales could help to understand why specific cities, or even neighborhoods within cities, are most susceptible to flood impacts.The mean-flow theory could also support the planning of new developments in anticipation of flood exposure, and the development of flood resilience plans and programs.Hence, the meanflow theory offers a promising entry point for urban planning of safer and more resilient cities.

Urban form data
We assimilate publicly available building footprints, digital elevation models (DEMs), and city boundary data sets for twenty cities worldwide, see Supplementary Notes I.For cities in the United States, building footprint data is obtained through Microsoft Maps' nationwide open building footprints data set 48 .For most cities outside the United States, we extract building footprints from OpenStreetMap 49 , a publicly available crowd-sourced database.To supplement this database, we also obtain building footprints for African and Southeast Asian countries through Open Buildings, a building footprint database developed through deep learning of satellite imagery 50 .We obtain the elevation data from MERIT 51 , a high-accuracy global digital elevation model at roughly 90-meter resolution.While the MERIT DEM was found to contain an RMSE of ~4m 52 , this error translates to a ± 0.008 change in slope and is likely present in all cities, having little effect on our overall findings.At a finer resolution of 30 meters, FABDEM 53 was also considered but was ultimately found to contain building heights that skewed the terrain elevations.We delineate the urban region based on the administrative boundaries of each city.

Synthetic urban forms
We conceptualize city urban form with rectangular landscape panels of length (L), width (W), slope (α) and Manning coefficient (n M ) populated with N square buildings of side length (D).Ensembles of urban landscapes representative of a specified porosity ϕ and order χ are generated through a hybrid reverse Monte Carlo (HRMC) algorithm, which was first introduced to improve upon the modeling of amorphous carbon structures 35,36 .Using user-defined constraints to inform the placement of buildings and an energy penalty to avoid unrealistic or physically improbable configurations, this probabilistic and iterative algorithm is capable of generating ensembles of urban forms with the targeted porosity and order (See Supplementary Note II and Supplementary Fig. S1 for additional detail).
Urban landscape panels were constructed with a length L=10,000 m, W=500 m, building size D=15 m, and a range of slopes, α=0.001-0.1.The channel was configured sufficiently long to balance gravitational effects and flow resistance from bottom shear and form drag, sufficiently wide to avoid edge effects, and sufficiently small to make fine-resolution modeling feasible over thousands of scenarios.Manning's roughness (n M ) is 0.02 s/m 1/3 in all simulations, representing a concrete surface, and an ensemble of M=3 realizations are utilized for each combination of porosity, order, and slope.We confirmed that a set of three simulations was sufficient for estimating the mean flood depth by checking sixteen different urban forms using M=30 realizations and seeing only a negligible difference.

Flow simulation
Flood flow for each realization of urban form is simulated using Par-BreZo, which uses a Godunov-type, finite volume scheme to solve the full two-dimensional shallow-water equations on an unstructured grid of triangular cells constrained by building walls 23 .Hence, buildings are represented by the so-called building-hole method 28 .
Every realization of urban form is independently meshed through constrained Delaunay triangulation using the Triangle mesh generation software 54 .To ensure a quality mesh, all cell angles are restricted to greater than 20 ∘ , and a spatially-dependent maximum area constraint on cell size is utilized.Within the upstream section of the urban form where we seek to achieve a uniform flow in which gravity and flow resistance balance each other, i.e., 0 < x < L/D = 33.3, the maximum area constraint is A c /D 2 = 0.02, where A c is the maximum area of a triangular cell.Hence, we use a cell size that is no larger than 2% of the size of a building.In the downstream portion of the urban form, i.e.,L/D = 33.3< x < L/D = 666.7, the mesh coarsens to a maximum cell area of A c /D 2 = 0.13 to reduce computational costs without impacting precision.See Supplementary Fig. S2 for a visual of the mesh.
To capture uniform flow within the upstream portion of each urban form, a flow rate of Q is evenly spread over the upstream edge of the domain, the downstream edge is treated as a free-overfall boundary condition, and the sides of the domain are treated as freeslip walls.The model is initialized with an initially dry floodplain and integrated into a steady state with an adaptive time step that maintains a Courant number of 0.8.Upon convergence, the uniform depth and velocity distribution are evaluated based on the flow conditions within the upper 100 m of the channel.For each model, the urban form and flow attributes ϕ, χ c n , Q, l c , and α, along with the resulting flow height h, velocity u, and their product hu, are provided in the spreadsheet 'Supplementary Data 1'.
We note that the domain length and width specification, grid refinement criteria, and convergence conditions were checked from an uncertainty perspective, and sensitivities were found to be negligible.Moreover, our multi-step simulation workflow, including (1) construction of synthetic urban forms, (2) mesh generation, (3) shallow water simulation, and (4) data post-processing, is implemented as a high-throughput framework.This framework runs on a highperformance computing cluster and streamlines the execution of parametric studies such as this one; see Supplementary Note III for details.

Global urban flood hazard analysis
Each of the 20 cities identified for this study is decomposed into 1 km x 1 km panels, Fig. 5a.While many cities have street drains with O(100m) spacing to collect rainwater, these systems are only designed to contain flows from relatively frequent events (e.g., ten-year return period) and thus substantially longer overland flow runs occur with the most severe events.The use of a value of 1 km allows us to resolve both the amount of overland flow and the slope of the ground, thus providing an assessment of urban flooding at the neighborhood scale.Moreover, we find that the trends in flood hazards revealed by the mean-flow theory do not depend on a specific value of the rainfall accumulation length scale.Scarcely developed regions with a porosity higher than 90%, and within city boundaries, are disregarded.For each panel, urban porosity and average building side length are calculated from the building footprint data.Additionally, the DEM of each panel is analyzed to find the maximum slope and direction.We subsequently calculate the effective chord length in the direction of maximum slope in each cell, Fig. 5b and Supplementary Table S1 provide the complete set of urban attributes computed for the twenty cities studied in this work.
The Emergency Events Database (EM-DAT) 45 is cross-referenced for "flash flood events" and "tropical cyclones" since 2000, which lead to direct runoff from impervious surfaces into the urban form and the concentration of flooding.A total of sixteen flash flood events are identified, and five of these contain estimates of monetary damages, see Supplementary Table S2.Additionally, a total of seven tropical cyclone events are identified and all contain information about monetary damages.For each event, we collect precipitation data from the PERSIANN system 55 .Precipitation data is reported hourly for the extent of each event, and the highest rate is applied to the city area to produce flow rate Q using the rational method and assuming that all precipitation becomes runoff (i.e., no interception or infiltration).
We do not consider flooding events that occurred before the year 2000 because PERSIANN satellite data for systematically estimating precipitation is not available, and because of substantial urban growth over the past twenty years.We also note that monetary losses recorded in the EM-DAT database are usually not confined to a city's boundaries but include all regions affected by the event.This may cause an over-or under-estimation of damages reported for a city.We also note that monetary losses from tropical cycles may be caused by multiple hazard drivers including wind, storm surge, and rainfall, whereas monetary losses from flash floods are mainly the result of intense rainfall 2 .

Fig. 1 |
Fig. 1 | Urban form across the cities in North and South America, Europe, Africa, and Asia is highly varied.Selected areas from (a) Boston, (b) Virginia Beach, (c) Chicago, (d) London, (e) Jakarta, (f) Sao Paulo, (g) Tokyo, and (h) Lagos.The building footprints, streets, highways, parks, and rivers are displayed in gray, yellow, orange, green, and blue, respectively.The coordinates of each city are defined along the top and left sides of the panel.Basemaps provided by OpenStreetMap (https://www.openstreetmap.org/copyright)used in conjunction with Microsoft's Building Footprints and Google's Open Buildings, data made available via an Open Database License (https://opendatacommons.org/licenses/odbl/).

Fig. 2 |
Fig. 2 | The qualitative study of flow pathways in synthetic urban forms shows the impact of urban porosity, order, and chord length on flood velocity distribution.Non-dimensionalized flood velocity maps for unit cells of urban forms with (a-c) square and (d-f) hexagon-like patterns with porosity ϕ=0.6 and Mermin order χ=1.0, 0.8, and 0.6, respectively.The color bar on the right represents the velocity scale for all velocity maps.Gray squares represent flow obstructions.While the unit cell l u is ~33D, the channel length is set to 20l u to mitigate size effects.The normalized chord length distribution Pðl 0 c Þ for each urban form is placed below each velocity map in panels (g-l); see the text for the definition.

Fig. 3 |
Fig. 3 | Considering mass and momentum conservation laws to develop a nondimensionalized relation between flood height, driving forces, and urban form attributes. a A schematic of the control volume through which water, shown in blue, flows between gray square obstacles.In dense urban environments, the floor shear is negligible compared to the other forces, and the gravitational force balances the form drag at the steady state.A limited number of one-factor-at-a-time flow simulations corroborate with the proposed dimensionless relationship, showing the flood height (b) scales linearly with flow rate per width (q), and (c) is proportional to the inverse square root of the bottom slope (α −1/2 ).

Fig. 4 |
Fig. 4 | Dimensional analysis collapses simulation data across urban forms (0.4 ≤ ϕ ≤ 1, 0:3 χ c n 1, 0.001 ≤ α ≤ 0.1) into a dimensionless flood depth, h -hi ffiffiffiffiffiffi ffi gD p =q.This depth depends on (a) separate functions of slope, porosity, and order for rectangular and hexagonal symmetry and (b) a single function of slope, porosity, and average chord length.

Fig. 5 |
Fig. 5 | The application of mean-flow theory to probe flood hazard globally.Flood hazards across cities worldwide are examined by (a) Subdividing cities into 1 km x 1 km tiles (shown in Vancouver, BC, Canada) and evaluating the topographic slope α, porosity ϕ, and order χ and (b) estimating the average chord length in the direction of descent to reveal (c) the porosity and chord length across 20 twenty cities; error bars show the margin of error with 90 percent confidence.Application of the mean-flow theory reveals: (d) the flow rate-normalized flood inundation H * = h hi=h ref and non-dimensionalized flood intensity P * = h hui=hu ref for twenty cities with constant flow rate of Q = 1m 3 /s, (e) the theoretical flow rate-normalized flood inundation H * and non-dimensionalized flood intensity P * for recorded flood events.Finally, (f) The normalized monetary damages L are compared to the product of the inflow rate q and urban form factor F ðϕ, l c , α, DÞ among different flash flooding events.In panels (c-f), the markers' colors represent the average bottom slope of each city defined by the color bar to the right side.Building footprint data for Karachi and Mumbai may not be sufficient.Basemaps provided by OpenStreetMap (https://www.openstreetmap.org/copyright)used in conjunction with Microsoft's Building Footprints and Google's Open Buildings, data made available via an Open Database License (https://opendatacommons.org/licenses/ odbl/).